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The nature of the low energy band of the Fenn a- Matt hews- 01 son complex: 
vibronic signatures 

Felipe Caycedo-Soler, Alex W. Chin, Javier Almeida, Susana F. Huelga and Martin B. Plenio 

Institute of Theoretical Physics Albert- Einstein- Alice 11 D - 89069 Ulm, Germany 

(Dated: 22 May 2012) 

Based entirely upon actual experimental observations on electron-phonon coupling, we develop a theoretical 
framework to show that the lowest energy band of the Fenna-Matthews-Olson (FMO) complex exhibits 
observable features due to the quantum nature of the vibrational manifolds present in its chromophores. 
The study of linear spectra provides us with the basis to understand the dynamical features arising from 
the vibronic structure in non-linear spectra in a progressive fashion, starting from a microscopic model to 
finally performing an inhomogeneous average. We show that the discreteness of the vibronic structure can 
be witnessed by probing the diagonal peaks of the non-linear spectra by means of a relative phase shift in 
the waiting time resolved signal. Moreover, we demonstrate that the photon-echo and non-rephasing paths 
are sensitive to different harmonics in the vibrational manifold when static disorder is taken into account. 
Supported by analytical and numerical calculations, we show that non-diagonal resonances in the 2D spectra in 
the waiting time, further capture the discreteness of vibrations through a modulation of the amplitude without 
any effect in the signal intrinsic frequency. This fact generates a signal that is highly sensitive to correlations 
in the static disorder of the excitonic energy albeit protected against dephasing due to inhomogeneities of the 
vibrational ensemble. 

PACS numbers: 82.53.Ps;87.14.E-;87.15.A-,87.15.H-;87.15.M-;87.14.E- 



I. INTRODUCTION 

The dynamics of photosynthetic complexes, such as the 
Fenna-Matthews-Olson (FMO) complex, have received 
rapidly increasing attention as an example of how nature 
appears to favor an intermediate noise regime, as sug- 
gested by the discovery of experimental evidence for long- 
lasting coherent oscillations present in signals arising 
from non-linear spectroscopy protocols^, and the the- 
oretical discovery that environmental fluctuations play 
a crucial role in explaining the very high efficiencies for 
excitation energy transport in such complexes^ 6 . The 
FMO complex is a trimer, with individual subunits com- 
posed of eight strongly coupled chromophores that serve 
as a link between the Chlorosome antennae and the re- 
action center of green sulfur bacteria^. Chromophores 
have a spatial extension that supports the existence of 
intramolecular vibrations^—. These vibrations play an 
important role when the chromophore structures are per- 
turbed due to the electronic charge redistribution upon 
photon excitation. 

Non-linear spectroscopy has proved to be useful in or- 
der to unveil the dynamics involved in excitonic trans- 
fer of light harvesting complexes, due to the fact that 
it is sensitive to excitonic quantum superpositions, i.e., 
excitonic coherences. Two dimensional (2D) non-linear 
spectroscopy can resolve the third-order polarization of 
the electronic system, from the heterodyne detected spec- 
trally resolved signal S(t\, £2, W3), arising from photo- 
excitation of three consecutive pulses with wave vectors 
ki , k2 and k3 , separated by time intervals t\ and t<^. 
Experimentally, it is customary to Fourier transform the 
time t\ dimension to yield S(ui, £2, W3) in order to gen- 
erate a 2D spectra parametrized by the waiting time t^. 



This spectra allows inquiring into the dynamics of ex- 
citation transfer in photosynthetic light harvesting com- 
plexes. 

The appreciation of the role of discrete vibrational 
modes on the excitation transfer process has evolved over 
timei^, and additional signatures of their presence in the 
2D signal will contribute to a thorough understanding. 
The main purpose of this paper is to show that the exper- 
imental evidence points towards the presence of a discrete 
vibronic structure in the lowest energy band of FMO, and 
that the dynamical contribution of these vibronic modes 
can be extracted from the information available through 
non-linear spectra protocols. 

Particularly, the multiplicity of states resolved in the 
low energy band in the FMO has been the subject of 
intense debate. It has not been possible to date, to 
uniquely assign its nature to electronic contributions. In 
this manuscript we propose a model for the lowest en- 
ergy band in the FMO that includes explicitly the most 
significant resonances of its environment spectral den- 
sity, which is consistent with the available experimen- 
tal evidence achieved through hole burnin g 11 ' 12 ! 15 , flu- 
orescence line-narrowin g 12 ' 16 , and accumulated photon 
echoil. We explore the possibility to detect the cou- 
pling of intra-molecular low energy vibrations of chro- 
mophores, to exciton dynamics, and shed light on the 
nature of the lowest energy band of the FMO (to be re- 
ferred as B825 from its absorbance maximum). We show 
that the experimental evidence can be explained with the 
existence of a representative vibronic resonance with en- 
ergy ~ 36 cm -1 and, beyond current suppositions that 
neglect higher order contributions, its overtone. Interest- 
ingly, this resonance has also been specifically addressed 
in resonant Raman spectra of the special pair BChla in 
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purple bacteria reaction centers^, inferred through the 
Stokes shift in LHC-II protein complexes Chla pigment 
of higher plants 19 , and exists within the energy range ex- 
pected for low-frequency vibrations of helical structures 
in protein o 20 i 21 . This strongly suggests that its origin 
are either the absorbing chromophores or a ubiquitous 
protein feature such as the a-helices both of which, are 
therefore not specific to the FMO structure. 

In order to introduce our model we comment on the 
outcome of different experimental examinations in the 
B825, and qualitatively reproduce their results with the 
assumption of the presence of a single excitonic state, a 
discrete vibronic dynamical contribution, and its over- 
tone. The dynamical contribution of this vibronic reso- 
nance is studied through the phase difference among sig- 
nals obtained in the neighborhood of an electronic diag- 
onal peak in the 2D spectra, which allows to distinguish 
the fundamental from higher vibrational harmonics. By 
unraveling the underlying vibrational dynamics in the 
2D spectroscopy time domain, we study the characteris- 
tics of the beating signal that arise due to this electron- 
phonon interaction (the phonon term will be used for 
the intramolecular vibrations of our interest). Accord- 
ingly, we explore the dynamics of non-diagonal contri- 
butions in the 2D spectra, to show that a modulation 
of the amplitude in the excitonic coherence signal arises 
due to these low frequency vibrational modes. We ob- 
tain an analytical expression for the excitonic coherence 
that explains this beating pattern as a result of vibra- 
tional wave-packet motion. We highlight that the co- 
herent oscillations are extremely sensitive to the degree 
of correlation in the static disorder of neighboring chro- 
mophores and propose a physical picture involving the 
spatially extended a-helices of the protein scaffold, that 
is able to explain the observed long lifetime for excitonic 
coherence. On the other hand, the modulation of the 
signal amplitude due to vibrational wave packet motion 
does not alter the intrinsic frequency of oscillations of the 
polarization signal in 2D spectra and therefore, produces 
no dephasing due to inhomogeneities of the vibrational 
manifold. The latter effect allows us to explain recent 
experimental results where, even though structural mod- 
ification in the FMO chromophores were accomplished, 
the dephasing rate of the electronic coherence in such a 
complex remained unaltered 2 ^. 



A. Brief history of the FMO low energy band B825 

Due to the energetic difference between the first exci- 
ton and higher lying excitonic states, the lowest energy 
band in the FMO complex has become the optical tran- 
sition in which to study and understand, the extremely 
long lifetime of excitonic coherences found in this pho- 
tosynthetic complex- 1 -. Nevertheless, the results obtained 
with different experimental techniques have not reached 
a consensus on the nature of this energy band. 

The first hole burning and fluorescence line-narrowing 



experiments performed in FMO 12 , showed two compo- 
nents in the B825, whose absorption differed by Alj m 
30 cm -1 . The possibility that these frequencies were 
zero-phonon and phonon side bands was considered but 
ruled out at the time, because the energy difference be- 
tween emission and absorption maxima in the B825 was 
AE ~40 cm -1 , and could not be reconciled with the 
calculated Huang-Rhys factor s = 0.3 from measured 
absorbance difference, that resulted in a Stokes shift 
2su = 2(0.3)30 fa 20 cnT 1 ^ AE. Rather, these lev- 
els were regarded as excitonic states arising from delo- 
calized excitonic components of the FMO trimer- 2 , a na- 
ture corroborated by theoretical calculations of the whole 
FMO trimer™ that suggested the existence of equivalent 
BChl's 7 (using the numbering scheme of Refii) on differ- 
ent units as the origin of the three excitonic states on the 
B825. Even though linear dichroism confirmed dereal- 
ization on B825 due to the broad linewidth in triplet- 
minus-singlet spectra 24 -, the linear dichroic absorption 
with respect to a continuous microwave field showed that 
the angle among the triplet and Q y transition dipoles in 
the B825 is constant throughout the band 2 ^, and close to 
that observed in single BChlc 26 . Since the triplet and sin- 
glet orientations are sensitive to excitonic derealization, 
the aforementioned equivalence among the B825 and sin- 
gle BChla, suggested that such a band was composed of 
localized excitations within a single BChla. This last 
turn of events was even further corroborated when the 
simultaneous fitting of absorption spectra, triplet-minus- 
singlet and linear dichroic absorbance, concluded that the 
B825 was composed of localized excitons on the BChl 3^1. 
Shortly afterwards, Stark hole burning studies^? com- 
pared the purple bacteria B800 band (composed of al- 
most completely localized excitonic states) with the B825 
to provide further support that in the B825, the inho- 
mogeneous broadening among equivalent BChls over the 
trimer generates excitons localized on single monomers 
of the FMO. The low energy band was successfully fitted 
with three gaussians, that were thought to represent the 
contributions of the localized lowest exciton states of each 
subunit of the FMO trimer^. The view that the contri- 
butions in B825 would come from energetic inequivalence 
arising from structural heterogeneity was also supported 
by the absence of excitonic dipole moment redistribution, 
probed by polarized hole burning for this energy band 11 . 
In that work, two important results were exposed. First, 
a greater homogeneous linewidth was found with greater 
energies within the B825, and confirmed the trend found 
by accumulated photon echo, of higher dephasing rates as 
the involved pulses were blue-shifted—^. Second, a very 
rich structure of satellite low energy holes regarded as 
intramolecular vibrations (with the biggest components 
having energies of 36 and 72 cm -1 ) in the B825 was first 
observed. 

The trend of increased linewidth at higher energy 
in the B825 was attributed then, to cascade excitation 
Forster energy transfer among individual chromophores 
on different subunits. However, Monte Carlo simulations 
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of the Hamiltonian under the hypothesis of these three 
different chromophores being BChl 3^2, were not able to 
simultaneously fit the lineshape and the gradual varia- 
tion of the homogeneous linewidth and dephasing times 
observed. Almost concurrently with hole-burning ex- 
periments, the technique of fluorescence line-narrowing 16 
successfully determined the phonon spectral density and 
confirmed the hole-burning experiments^ about the ex- 
istence of discrete phonons with frequencies 36 and 70 
cm -1 moderately interacting with excitons. The possi- 
bility that multiple dynamical components on this band 
could arise from these modes, was first suggested. 

Difference fluorescence line narrowing^ provided fur- 
ther evidence of three contributions on the B825, and 
suggested once again, the possibility of exciton der- 
ealization in the lowest exciton state since new simula- 
tions of the B800 in purple bacteria harvesting structure, 
showed a derealization length across 2-3 pigments 31 . Re- 
cently, molecular dynamics simulations^ showed an ap- 
preciable narrowing of the of site energies in a trimer 
compared to a monomer. This narrowing renders infea- 
sible that static structural inhomogeneities in the trimer 
are the cause of the B825 multiple dynamic contributions. 

At present, the consensus view is that of a sin- 
gle excitonic state in the low energy band to describe 
the FMO optical spectra. This view was corroborated 
with more detailed calculations of the FMO electronic 
Hamiltonian 3 ^, that have shown increasing success in 
fitting optical spectra, when linewidth assignments of 
such transitions involved a fit of the spectral density 
coarse grained envelope with a global Huang-Rhys fac- 
tor s ~ 0.5. However, no model has been able explain si- 
multaneously the facts that within the B825, 1) there is a 
single electronic excitation where at least two representa- 
tive resonances appear, 2) a marked trend of greater ho- 
mogeneous linewidths (or dephasing rates) while probing 
higher energies was found. Beyond these points (high- 
lighted in a recent review^), the fact that 3) a very 
rich intramolecular vibrational structure is observed but 
never explicitly accounted for. This suggests the lack of a 
complete understanding of the low energy band of FMO. 



II. THE MODEL 

Our immediate purpose is to explain the possible na- 
ture of the multiple resonances found with increasing 
linewidth for blue-shifted probing of the B825. At the 
same time, the model must be consistent with the ap- 
proximate 30 cm -1 absorption energy difference among 
these transitions, with the observed Stokes shift and 
with the fluorescent line narrowing spectral density. The 
dipole moment of such transitions must not exhibit vari- 
ations in polarization while scanning this energy band. 
Let us examine these issues in detail. 
The most recent Huang-Rhys factor s — 0.5 describing 
the global electron- vibration interaction, is better able to 
explain the Stoke's shift of 40 cm -1 observed, occurring 



mainly due to the 36 cm' 1 mode, i.e., 2sw = 2(0.5)36 ~ 
40 cm -1 . The whole phonon wing Huang- Rhys factor is 
an upper bound for the interaction strength among this 
mode and electronic excitations. In particular, hole burn- 
ing spectral assigned Huang- Rhys factor s = 0.12 for 
the most representative resonances, namely, resonances 
having energies of 36 and 72 cm -1 . Fluorescence line 
narrowing 16 reported for both these resonances, a Huang- 
Rhys of s = 0.01, on a procedure that might overestimate 
the weight of the phonon wing. As the authors explic- 
itly point out, it is heuristically assumed that individual 
modes just refine details over the phonon wing, and there- 
fore minimally contribute to build up the maximum of 
the spectral density function. The result of either model 
differ beyond the purely practical aspect of a fitting pro- 
cedure, since a discrete contribution may lead to coherent 
evolution that will result in crucial features to be studied 
below. Therefore, if the coupling of this mode is up- 
per bounded by the global Huang-Rhys factor s = 0.5, 
and lower bounded by the reported hole-burning Huang- 
Rhys value s — 0.12, a perturbative analysis is dubious. 
Moreover, the vibronic contribution does not represent 
a redistribution of dipole moment, allowing a single po- 
larization throughout the B825. Hence, we proceed with 
the explicit inclusion of this mode. 

The Hamiltonian for electron-vibration interaction 
reads 

H ° = E ( t + ^~ k ^ hk + b+ ^) a * + ^ h w 

k 

where &i and are usual spin 1/2 Pauli matrices, and 
annihilation (creation) operators of the fcth boson field, 
the electron- vibration coupling strength y/stcuk is given 
in terms of the phonon frequency u>k and Huang-Rhys 
factor Sk- 
in the present section, our aim is to consistently de- 
scribe the homogeneous broadening trend as the prob- 
ing laser field is blue shifted when vibrations are explic- 
itly included. To that end we consider a single vibration 
and the laser interaction with the charge distribution of 
the chromophore under the rotating wave approximation 
(RWA) in a frame rotating at the laser frequency, we 
obtain a hamiltonian: 

H = H + (A - uj a )a z + ^a x (2) 

where A = luq — lol and f2 are laser detuning from the 
energy of the zero phonon line (ZPL) electronic transition 
and Rabi frequency, respectively. 

Now, the broadening of the resonances is of utmost im- 
portance in the historical discussion and is required for 
the analysis of spectra, as it influences the interaction 
of the system with the environment. The homogeneous 
linewidth of the resonances will primarily depend on the 
different processes inducing decoherence beyond the en- 
semble average. Accumulated photon echc-^ fitted the 
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dephasing rate in 1.8-50 K range with the expression: 
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(3) 



for the homogeneous dephasing time T2, and highlight 
the main processes involved: population relaxation, pure 
dephasing from spectral diffusion promoted by the sur- 
rounding protein, and a thermally induced activation 
process. Given that the FMO complexes used in optical 
experiments lack the RC, the first term in eq.([3]) is lower 
bounded by the fluorescent relaxation. The second term, 
is attributed to dynamics of tunneling among energy min- 
ima of different structural configurations in the protein, 
and characterized by a rate that obeys a power law tem- 
perature dependence oc T 13 , typical of glassy hosts. The 
third terms involves thermally activated transitions not 
specified ini£, that accounts here for the thermal equi- 
libration of the vibrational manifold. Note that eq.© 
only describes thermal relaxation, and therefore is best 
suited to the red edge of the B825. A fit to the red 
edge photon echo signal was accomplished with free pa- 
rameters {T 2 ,j Z7 j,SE} = {420 ps,0.09 ns" 1 /^ 1 ' 3 , 23.6 
ns -1 , 15 cm -1 }. In our model the most prominent ther- 
mally activated process occurs among states whose en- 
ergies involve boson quanta of energy 36 cm -1 which 
restricts SE = 36cm" 1 . A good agreement among 
both sets of parameters is accomplished for a constant 
7 = 70.09ns~ 1 =2.36 cm -1 (results not shown). 

In a first model, these mechanisms can be described 
under the Born-Markov approximation with super- 
operators d(p) = LipLf - j(L+Lip + pL+Li) in a 
Lindblad-type master equation (h = 1): 



d tP =-i[H,p] + J2^(p) 



(4) 



with L\ = vTVr for population relaxation (a ^ are 
lowering (rising) Pauli spin ladder operators) and L2 = 
\/2{n + 1)76, L3 = ^/2n r yb + describing thermalization 
to temperature T = 1/fik of the mode having an equi- 
librium average quanta n = (er^ — 1) _1 — . The factor 2 
in this latter pair is introduced since the dephasing pro- 
duced is half the rate of relaxation/absorption processes. 
The protein configuration tunneling induced dephasing is 
accounted by the operator L4 = y / j^T^a z . The magni- 
tude of rates T, 7 and j z are hence provided from the de- 
phasing rate obtained through accumulated photon echo 
technique^. 



III. RESULTS 

A. Linear spectra 

The fluorescence intensity in the line narrowing sig- 
nal is proportional to the excited state population in 
the stationary state of the chromophore (cr + cr _ ) t _ i . 00 , ob- 
tained from eq.(TJ|, and presented in FigQJa). First, note 



that with both the upper and the lower bounds of the 
Huang-Rhys factor, three fluorescent transitions (high- 
lighted by arrows) can be resolved: the zero-phonon line, 
the phonon side band, and the overtone, which coincides 
with the 70-72 cm -1 resonance reported inii. It should 
be noted that there is no need to include another vibra- 
tion 70-72 cm -1 since a single mode and its harmonic 
lead to the most representative resonances at ~ 36 and 
~ 72 cm -1 . Increasing the Huang- Rhys factor to s = 0.5 
results in a higher overlap between levels having differ- 
ent boson quanta in the excited and ground electronic 
states, and therefore a redistribution of fluorescence in- 
tensity increases the contribution of higher harmonics. 
It should be noted that even a moderate Huang-Rhys 
factor s = 0.12 results in three dynamical contributions 
at T —77 K, as shown in the inset of Figflja). At this 
temperature, the homogeneous phonon side band is com- 
mensurate with the zero-phonon line, and therefore, even 
in the presence of ensemble inhomogeneities, both tran- 
sitions should be of importance in the description of the 
FMO in 2D spectroscopy (T =77 K, Ref-i). 

Having captured these three dynamical contributions, 
we proceed with the analysis of the individual linewidths. 
The hole burning experiments, sample the absorbance 
A(cj) proportional to the Fourier transform (FT) of a 
two-time correlation function, explicitly: 



A{uj) oc lim 
t— > 00 



(<J-{t)a+(t + r))e^^ r dr, (5) 



which is presented in Figfljb) (calculated with the quan- 
tum regression theorem according to the Lindblad form 
master equation eq.(|4])) for the relevant burning frequen- 
cies. In this figure, it is shown that higher energy on 
the burning beam generates greater linewidths of the 
holes in the absorption spectrum. With this model, the 
linewidths, say 70, 71, 72, of the ZPL, the phonon side- 
band, and overtone, respectively, show a difference that 
persists for the whole range of temperatures measured 
as shown in the inset of this figure. This is in com- 
plete agreement with photon echo experiments (1.8-50 
Kfil). The broadening can be understood from transi- 
tions involving n-th vibrational states, that become more 
populated with increasing temperature and whose ther- 
malization rate is proportional to the excitation number 
(£2(3) \ n ) ( n \) n - The model also captures two impor- 
tant facts found in hole burning and photon echo experi- 
ments. Firstly, the greatest difference between linewidths 
occurs at low temperatures with a trend to converge as 
temperature is increasedAi. Secondly, before convergence 
is reached to equal linewidths, the three components are 
widening in the same fashion 17 . This behavior illustrates 
the prevalence of different mechanisms depending on the 
temperature. At low T, spontaneous emission and ther- 
malization decay are the key mechanisms, the former be- 
ing independent of temperature. At intermediate tem- 
peratures, dephasing induced by protein conformations is 
equal in any of the transitions and represents the trend of 
the linewidths to equalize. At high T, the onset of ther- 
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ACcm" 1 ) a>-cji[cm ] 

FIG. 1. In (a) fluorescence intensity from solution of stationary state in eq.Q. Dashed and continuous correpond to Huang- 
Rhys factors s — 0.5 and s = 0.12, at T=4K, drawn with the same scale. The arrows point the frequency of the zero-phonon (0 
cm" 1 ), phono side-band (36 cm -1 ), and overtone (72 cm -1 ). The inset shows the fluorescence intensity at T=77K, s = 0.12. 
In (b) A{cj) is normalized for the three resonances and correspond in continuous, dashed and dotted, to the A = 0, 36 and 72 
cm -1 resonances, that have a FWHM of {70,71,72} = {1.6,3.5,5.6} cm -1 at 4 K, from a fit with lorentzian functions. Inset 
(b), shows the variation of the lorentzian widths with temperature for the ZPL (70, circles), sideband at 36 (71, boxes) and 
excited harmonic 72 cm" (72, diamonds), that yields to values of {70, 71,72} = {6.2,8.1, 11.2}cm" 1 at 77 K. If otherwise not 
stated T = 7.9 x 10~ 2 cm -1 , 7 = 2.36 cm -1 , -y z0 = 3 x 10" 3 cm" 1 , s = 0.12 , u = 36 cm" 1 , Q = 10" 3 cm -1 . 



malization as the dominating dephasing mechanism in- 
duces equal increment in all linewidths with temperature. 
In practice, the reddest photon echo signal according to 
eq.([3]) from accumulated photon data^i, extrapolates to 
a dephasing rate of 3.3 cm" 1 . Based in our theoretical 
results, this rate is « 7o/2 = 6.2/2 cm -1 (see the inset 
Figfrjb) at T=77 K). The difference of the data with our 
ZPL linewidth can be traced back to mixing of echoes 
originating from the ZPL and the side-band within the 
ensemble. 

At this stage, we conclude that by inclusion of a vibra- 
tional mode with frequency of 36 cm -1 and an Huang- 
Rhys factor s > 0.12, we are able to capture three dy- 
namical contributions with equal polarization, whose po- 
sition in absorption spectra and linewidths trend follow 
the experimental observation in the lowest energy band 
of the FMO. Without inclusion of a such mode, the ab- 
sorptive response showing three dynamical contributions 
below T=10 K 11 cannot be explained relying solely on an 
inhomogeneously broadened single excitation electronic 
contribution. 

Additional information is available in the 2D non- 
linear spectra that might better characterize this low en- 
ergy vibration, and further support (or reject) its role in 
the additional dynamical component on the lowest en- 
ergy band of the FMO. 



B. Non-linear 2D spectroscopy 

1. Background 

In 2D spectroscopy, the heterodyne detected signal 
S(u)i,t2,wa) allows straightforward comparison with the 



FT of the theoretical analogue 

5(wi , ta, w 3 ) = FT[2 Im Ri *a, h^J ] (wi , *2, w 3 ) 

(6) 

for a two level electronical system with energy difference 
(jj eg among excited and ground states |e) , \g) respectively. 
In this case, the response functions R2 and i? 3 (Ri and 
-R4) usually termed as rephasing (non-rephasing) paths, 
allow calculation of signal in the direction of observation 
k s = ki + k 2 + k 3 (k s = ki — k 2 + k 3 ) surviving the 
rotating wave approximation. In order to present the 
data in a single quadrant, the FT for the non-rephasing 
and rephasing signals have conjugate variables {— u>i, cj 3 } 
and {wi,u; 3 }, respectively. In detail, 

Rn = | M 4 | e ±i ^( tl T*3) e /..(t3,t2,ti) ( 7 ) 

where upper (lower) signs correspond to rephasing (non- 
rephasing) terms, respectively. In the above expression, 
/i is the electronic transition dipole moment and the func- 
tion f n (t^,t2,ti) depends on the lineshape function g(t) 
evaluated in combinations of these three time intervals^ 3 -. 
The lineshape function of a Brownian oscillator 

g(t) — isujt + s(i sin cut + coth(/3ftw/2)[l — coswi]) (8) 

describes the contribution to the response function from 
a discrete vibration with frequency a.' 13 . In this particu- 
lar case, the function g{t) generates through /„(i)2£, an 
increase of the excitonic frequency by the reorganization 
energy suj, uj eg — > uj eg + slo, produces a global reduc- 
tion of the signal by the Debye- Waller factor squared, 
e -2scoth(fScu/2) ^ anc j nag additional time dependent terms 
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of the form e s(i*™("t)-cath(0u/2)co S (ujt)) ^ The Fourier 
transforms required for the 2D plots are performed here 
analytically from a Taylor series expansion e° s ~ l+Os+ 
(Os) 2 + • • • . The exponential function has a rapid conver- 
gence on this expansion for s coth(/3w/2) <C 1, partially 
fulfilled for the Huang- Rhys factors of interest at the ac- 
tual temperature of experiments (0.12 < s < 0.5, lo = 36 
cm- 1 , T=77 K; scoth(/3w/2) = 0.37- 1.53). A first or- 
der expansion of this kind was developed in Ref<2&. In 
the following we further develop the expansion to second 
order, to study the possibility of resolving the higher har- 
monics (conspicuous in the absorption and fluorescence 
intensity) and develop a set of features able to better 
characterize vibronic resonances in the 2D spectra. 

The decay of the time resolved heterodyne signal in 
the impulsive limit, is equivalent to the linewidths of 
spectrally resolved observation under continuous- wave il- 
lumination. Henceforth, the effects of the bath in the 
2D spectra will be introduced in the aforementioned ex- 
pansion for terms proportional to e ±lu>tl , e ±lwt3 , e ±2lUJtl , 
e ±2 *"* 3 , with the Fourier transform argument +i"f\ 
on the former two, and lou$\ + i"/2 on the latter pair, to 
include the homogeneous broadening for side-bands and 
overtone. The terms lacking of any of these factors, are 
Fourier transformed with an argument + ijo mod- 
eling the width of the zero phonon line. The addition of 
these imaginary contributions in the Fourier transform 
argument is consistent with a decay to the stationary 
state in the waiting time domain due to the substitution 
giwta giwt 2 g— 71*2 and e 2i "' 2 — > e^^^e -72 ' 2 



2. Vibronic features in 2D spectra: Diagonal peaks 

Genuine vibronic sidebands in 2D spectra are chal- 
lenging to discriminate from the electronic contributions. 
Fortunately the exciton 1, in the FMO, has an en- 
ergy separation > 150cm -1 from all other excitonic states 
which allows its spectral resolution from the other elec- 
tronic contributions^ 3 -?. The 2D spectra can be used in 
order to characterize the most salient features of the vi- 
bronic spectra in this lowest energy band, such as the 
resonances found at 36 and 72 cm . 

The result of the second order expansion of the ho- 
mogeneously broadened 2D spectra for exciton 1 in- 
teracting with a vibrational mode at different waiting 
times t2, is presented in Figj2][a)-(b) for the rephasing 
and non-rephasing contributions, respectively. A hy- 
perbolic sine scale is used to highlight small features. 
First, note that in these figures, an asymmetry is promi- 
nent in the absorption dimension. This is highlighted 
in the zero-phonon emission profile (continuous line), 
where the difference in height of the peaks at W3 = io eg 
and cji = ui eg ± lo accounts for Stokes and anti-Stokes 
contributions due to differences among populations of 
the initial canonical state in the vibrational manifolds, 
p = J2 n e-P™\n)(n\/Tr{.}. 

Second, notice that the rephasing and non-rephasing 



contributions to the response function are sensitive to 
different processes on the vibrational manifold at a given 
waiting time. In detail, at zero waiting time the vibra- 
tional structure around the ZPL peak is prominent in the 
rephasing spectra, where even the overtone is conspicu- 
ous at lo\ = Co eg + 2cj, while such structures are almost 
inexistent in the non-rephasing signal at the same wait- 
ing time. For a waiting time ti = Tr/^j corresponding 
to half a vibrational period, the behaviors are exchanged 
leading to a non-rephasing contribution with very well 
resolved side-bands 3 -^ 

A better contrast of the side-bands is possible in homo- 
geneously broadened spectra, when instead of following 
the established procedure eq.® of adding the rephas- 
ing and non-rephasing responses (Figl^c)), we proceed 
with their subtraction (Fig[2][d)).The relative phases of 
rephasing and non-rephasing signals induce a stronger 
ZPL when added, but strongly suppress the weight of 
the ZPL to make the side-band structure gain relative 
weight when subtracted. This example allows to point 
out that the relative phases of the components that lead 
to the total signal are sensitive to algebraic operations 
over the individual contributions. However, this sce- 
nario where the average among the ensemble is not per- 
formed, neglects the static inhomogeneity of importance 
in the ensemble required to individually address rephas- 
ing and non-rephasing contributions in different spatial 
directions. 

In the following, the exciton energy inhomogeneities 
are studied, considering the second order expansion of 
the line width function. The localized nature of exci- 
ton 1 (at pigment 3) has been confirmed with numeri- 
cal simulations 39 , using a site energy distribution width 
of 110 cm -1 which due small exchange narrowing leads 
to a value of <je = 102 cm -1 for the total width of 
the B825 band. Recent numerical simulations of the 
FMO 2D spectra used a comparable width for the ex- 
citon 1 energy distribution (100 cm -1 ) 40 . The diagonal 
1-1 peak for rephasing and non-rephasing signals is pre- 
sented in Fig|3fa)-(b), respectively, where the inhomoge- 
neous broadening is introduced according to a gaussian 
distribution of the excitonic energy with standard de- 
viation cte — 102 cm . Interestingly, a very dramatic 
change occurs when moving from the homogeneous to the 
inhomogeneously broadened ensemble, where each path 
(rephasing/non-rephasing) is sensitive to addition of sig- 
nals arising from different processes. On the one hand, 
Fig|3Ja) shows the zero-phonon emission line (continu- 
ous) of the rephasing component, in order to highlight 
side-bands located at lo\ — lo$ = ±lo. On the other hand, 
the non-rephasing part, Fig[3jb), is sensitive to the over- 
tone, highlighted by the continuous lines over the surface 
at the ZPL and overtones at lo\ — L03 = ±2lo. A better 
contrast of the vibronic structure is accomplished in the 
non-rephasing 2D spectra due to the drastic reduction of 
the ZPL. Hence, the inclusion of inhomogeneous broad- 
ening does not simply regards to a general widening of 
the homogeneous resonances breath, since it affects in a 
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FIG. 2. In (a) and (b) we present the second order expansion for rephasing and non-rephasing paths, respectively, (c) the 
addition of rephasing and non-rephasing contributions, while (d) subtraction FT[i?2 + -R3 — (Ri + FU)]. In plots (a)-(b) an 
arcsinh scale has been used in order to discuss the side-band features, while (c)-(d) use a linear scale for purposes of comparison 
and t2 = 450 fs. In all plots the zero phonon emission line is shown in continuous line to guide the eye. Exciton 1 energy 
uj eg = Lo eg + £V SiUJi =12121 cm" 1 —, vibrational mode energy uj = 36 cm" 1 , and Huang- Rhys factor s — 0.12. 



different fashion to each of the resonances. In short, a 
discussion of this reduction will allow us to understand 
the qualitative differences among homogeneous and in- 
homogeneously broadened spectra. 

The inhomogeneous broadening in the excitonic fre- 
quencies corresponds to widening along the diagonal in 
the 2D spectra, with a ZPL maximum at ui\ — 0J3. We 

I 



have chosen the waiting time long enough in Figf3Ja)-(b), 
to avoid the ti dependence on the 2D spectra. These re- 
sults remain valid for different waiting times as will be 
demonstrated below. The expression for the zeroth or- 
der (waiting time independent) contribution of the ZPL 
reads: 



(7g + ("i-£W 2 ) a 



(9) 



where the upper and lower signs correspond to rephas- 
ing and non-rephasing expressions, respectively. The real 
part of the rephasing contribution is always positive for 
any value of the exciton 1 energy E % p 1 , while that of the 
non-rephasing part, is positive for E^ t — w\ < 70 and 
negative otherwise. The imaginary part of the rephas- 
ing signal is zero at the diagonal, while that of the non- 
rephasing, is positive for < uj\ and negative oth- 
erwise. This situation leads to have at a single point 
in the diagonal of the 2D spectra, only positive con- 
tributions for the rephasing part that will add up to 
build the ZPL peak. On the other hand, in the non- 
rephasing path the signals have varying signs depending 



I 

on the exciton energy, which when added up in the in- 
homogeneous sample, will lead to a decreased ZPL. This 
fact is quantified in FigJ3Jc)-(e), where the weighted con- 
tribution e~^ E *^~ Q "^ ^ 2(Te S(uji,lo 3i E^^ of the complete 
second order expansion as a function of the difference 
Efa —Cj eg are presented. The Figf3fc) shows that the sig- 
nal at LUi—uj3—uj e g has alternating (equal) signs for real 
and imaginary parts of the non-rephasing (rephasing) sig- 
nal. The trend continues (Figf3]Jd)) for points regarding 



the phonon side-band (o>i = uj eg + 



; ) where 



the real (imaginary) part of the rephasing signal presents 
over a great range of the whole inhomogeneous distribu- 
tion, positive (negative) values that constructively add to 
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FIG. 3. In (a) and (b), rephasing and non-rephasing 2D spectra of exciton 1 diagonal peak of FMO complex, broadened 
by excitonic energy inhomogeneities. In (c), (d) and (e) are presented the real (continuous, blue online) and imaginary 
(dashed, red online) contributions to the inhomogeneous signal at points (cji = CJ3 = cj eg ), (wi = u) eg + u),u)3 — oj eg ) and 
[u}\ = uj eg + 2u, tJ3 = uj eg ) respectively, for complexes whose energy differ E^ J1 — d) e9 from the ensemble average as explained in 
the text. In (f) and (g), rephasing and non-rephasing 2D spectra of exciton 1 diagonal peak of FMO complex, broadened by 
both excitonic and vibrational mode energy inhomogeneities. oe = 102 cm -1 , — 10cm _1 , — 2 ps. 



build up a appreciable side band in the rephasing path. 
However, rephasing and non-rephasing parts have posi- 
tive and negative contributions along the inhomogeneous 
distribution from the signals arising at the overtone reso- 
nance = a) eg +2aj, = uj eg ) which lead to a reduction 
commensurate in both. The concomitant reduction of the 
side band and drastic reduction of the ZPL for the non- 
rephasing direction, explains the surprising enhancement 
in the contrast for the signal arising from overtone, due 
to the (usually troublesome) sample inhomogeneities. 

The inhomogeneities are not only ascribed to the ex- 
citon energy, but also the vibrational manifold. For in- 
stance, hole burning can select a subensemble having an 
optical transition frequency, but the inhomogeneities in 
the satellite structure due to the vibrational manifold 
will not be different from those of any other subensem- 



ble. Therefore, the width of the satellite structure in 
these experiments reflects that of the ensemble. The low 
energy satellite structure in the B825 was observed in 
hole burningii, with a width of 10 cm -1 for the to = 36 
cm -1 resonance. Hence, we proceed to perform an aver- 
age of the vibrational frequency, according to a gaussian 
with cr w = 10 cm , and present the results for rephas- 
ing and non-rephasing paths in FigJ3{f)-(g), respectively. 
In the rephasing contribution, this average smears out 
the vibronic structure, while some of it remains in the 
non-rephasing path. Since the experiment involve macro- 
scopic samples, the vibrational manifold ensemble inho- 
mogeneities will lead to a reduced structure in the 2D 
spectra, and a simple inspection will not permit (in gen- 
eral) a straightforward inference of the presence of dis- 
crete vibronic contributions. Furthermore, the sampling 
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FIG. 4. In (a), Feynman diagrams involving the phonon- 
side band. Letters label electronic states, the superscript the 
boson state. In (b) are shown the waiting time domain sig- 
nals at the homogeneously broadened peak at oj\ = U3 eg + cj, 
cj3 = uj eg for the set {70, 71, 72} = {6.2, 8.1, 11.2} (left panel) 
and {70,7172} = {50,136,156} cm _1 (right panel). In (c) 
are presented the results for inhomogeneous broadened peaks 
with <te = 102 cm" 1 , a u = 10 cm -1 : in the left panel are 
shown the signals resolved in waiting time of the ZPL (top), 
side-band (middle) and second harmonic (bottom) peaks; in 
the right panel the slope <9 t2 Re[S']t 2= ioofe (long enough to 
avoid overlap with the second pulse ti > 16fs, the equality 
fulfilled for the pulses duration achieved in») is presented; 
arrows highlight the positions of the ZPL, side band and 
overtone. In all plots non-rephasing (continuous, black) and 
rephasing (dashed, blue) path signals. 



of the signal in t\ is done during a finite time that im- 
poses a minimum frequency, hence a maximum resolu- 
tion among frequencies. Currently, studies of 2D spectra 
tracking the vibrational wave-packet motion, used a scan- 
ning range for t\ from -100 to 100 fs, which allows a res- 
olution in Awi ~ 52 cm -1 for rephasing/non- rephasing 
signals, or half of it, Awi ~ 26cm" 1 for the total signal 38 . 
Accordingly, the discovery of vibronic structure in the 
ensemble becomes difficult relying on the inspection of a 
single 2D spectrum alone. However, the vibronic contri- 
butions can be detected from the waiting time £2, which 
induces features that allow us to infer the presence of 



discrete phonon modes. 

The additional structure whose peaks lie at uj\ — W3 = 
zhnu, imply rotations occurring in the vibrational mani- 
fold that, as will be seen, induce a phase shift in the signal 
on waiting time dimension. To illustrate our point, we 
start with the second order line-width expansion at the 
homogeneous level and focus on the real part of the sig- 
nal concerning the product of absorptive (non-dispersive) 
spectra in u)\ and LJ3. 

For instance, at w\ = Cj eg + lu, W3 = £b eg the signal is 
constructed from the Liouville paths contribution shown 
in the Feynman diagrams of Figf4ja) , where time evolves 
from bottom to top over the vertical lines that represent 
the ket and bra evolution. The wiggly arrows denote in- 
teractions with electromagnetic pulses which excite (de- 
excite) the system when they are directed into (out of) 
the ket-bra vertical lines. The additional kets or bras, 
are used to clarify the transitions involved after a given 
excitation. The diagrams R\ and R2 (R3 and R4) in- 
volve rotations in the excited (ground) state manifold 
for the waiting time domain. Observation of the point 
u>i = uj eg + lj is generated when the bra (ket) (e n+1 | 
(|e™ +1 )) is excited with the first pulse. On the other 
hand, at W3 = uj eg the evolution involves the same boson 
quanta in both bra and ket. Therefore a flip of the bo- 
son state n + 1 — > n is required in order to end up in a 
population ket-bra combination before the trace opera- 
tion is performed. For diagrams R3 and R4 this rotation 
can be accomplished when the second pulse de-excites 
the system to (g n+1 | and respectively. However, 

for i?i and R2 the rotation must happen in the wait- 
ing time, and therefore, the signal is sensitive to the the 
required flip in the vibrational manifold. This rotation 
can be explained along the mentioned mathematical pro- 
cedure, but can also be understood by considering the 
Frank-Condon principle. 

After pulsed excitation, the system is left in a super- 
position of eigenstates of the Hamiltonian (JT]) , which also 
represent a vibronic progression \h) whose energy differ 
by lo among consecutive levels in a given electronic man- 
ifold. As the thermal population involves the lowest en- 
ergy levels, the contribution to the signal will primarily 
arise from the ZPL and the phonon side-band, approxi- 
mately described by \<f>) ps a sin(wt) |5) + a\ cos(ujt) |l). 
Simply stated, the signal at u>U3> — Cj eg ± nu is mainly 
contributed by the \n) bosons state population, as a wit- 
ness of excitation and dynamics of the vibrational wave- 
packet motion. The phase of the signal will vary as 
the 2D spectrum is swept reflecting the phase difference 
among populations of the levels being scanned. In the 
waiting time, FigEJb) , the signal presents the discussed 
phase shift rs tt among rephasing and non-rephasing con- 
tributions (remind FigJ3fa)-(b)) . However, this panel 
also shows a phase shift among ZPL and side band sig- 
nals of 7T, to resemble the just mentioned population in- 
version. A population picture can be misleading for the 
phase difference among the side band and overtone, since 
this latter resonance develops a signal with strength pro- 
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portional to s 2 , small enough to be contaminated by the 
tails of the other resonances. The right panel of FigfJlb), 
shows that if the breath of the homogeneous broadened 
resonances is amplified in accordance with loosing the 
mode discreteness, the phase shift among different points 
in the 2D spectrum is lost when scanning the neighbor- 
hood of a diagonal peak. Hence, a smeared out interfer- 
ence of several indistinguishable contributions in a given 
area of the 2D spectrum reduces the phase shift differ- 
ences. 

In FigUJc) left panel, is shown the result of averag- 
ing static inhomogeneity to the signals resolved in the 
t-i dimension. Here it can be first noticed the reduc- 
tion of the ZPL and side band of the non-rephasing 
compared to rephasing path. Secondly, this procedure 
only slightly reduces the phase difference among rephas- 
ing and non-rephasing components. Third, the phase 
differences among resonances on a given path is per- 
sistent among ZPL, side band and overtone. A varia- 
tion in the phase leads also to a variation of the slope 
9t 2 Re[5](wi, W3) at a fixed t 2 , before the settlement of 
the stationary state. This variation is presented in the 
right panel of FigfJJc), and specially highlights (by fill- 
ing the area underneath) the change of sign of the signal 
slope for a fixed waiting time as a function of the coordi- 
nate wi . Even though the actual values of this phase are 
affected by the presence of the ensemble inhomogeneities, 
its variation remains, and since its origin is traced back 
to the discreteness of the vibrational manifold, it can be 
used as a tool to diagnose the presence of vibronic dis- 
crete structure beyond the unavoidable presence of static 



where the last two terms correspond to the electronic 
Hamiltonian, and Hu — Wjbf&j, Xi = bi + b\, bi (bf) 
are annihilation (creation) operators of the ith vibra- 
tion. Diagonalization of the electronic Hamiltonian leads 
to excitonic states \ipi) — -^=.(|3) — o|4)), ^2) = 

^=f(a|3) - |4)), parametrized by a = (V 8 2 + W 2 - 
5)/2V which for 6 = tu g 4 — w g 3 describes the limiting 
cases of localized (|V| -C S, a « 0) and fully delocalized 
(|V| S> 8, a s» ±1) excitons. 

In order to determine the effect of the vibrational modes 
in the exciton dynamics in the particular case of the FMO 
excitons 1 and 2, we proceed with a numerical solution 
constrained to the experimentally determined parame- 
ters. The dipole-dipole couplings in agreement with lin- 
ear absorption, circular dichroism and linear dichroism 33 , 
and the experimentally determined exciton transition 
energies 37 of the full FMO, were used in an iterative 
procedure as described in Refill to determine the site 
energies of 12142 cm -1 and 12315 cm -1 for sites 3 and 



disorder. 

The 2D spectra consist of diagonal and non-diagonal 
peaks. We have covered the behavior of signatures in 
both homogeneous and inhomogeneously broadened sam- 
ples for the diagonal peaks. Next, we consider a promi- 
nent signature that arises in the waiting time evolution 
of non-diagonal excitonic peaks. 

The single excitation peak on the model just presented 
is unable to capture the behavior of excitonic coher- 
ence. Moreover, it has been previously stateo-^Z, that 
extracting the dynamics of non-diagonal peaks in non- 
linear spectra, that is the excitonic coherences, has the 
advantage of less congested spectra. Therefore, it is an 
important contribution to understand the features that 
will arise in the cross peak dynamics on the FMO light 
harvesting complex. 



3. Vibronic features in 2D spectra: Non-diagonal peaks 

In the following, we investigate the excitonic coherence 
of a dimer interacting with localized vibrations in order 
to address the FMO exciton 1-2 cross peak, when the 
coupling to the lu = 36 cm -1 mode is accounted for. The 
sites with the greatest contribution on excitons 1 and 
2, are chromophores historically- 7 " named as sites 3 and 
4^ (states |3) and |4)). The Hamiltonian restricted to a 
single excitation for these two sites, each interacting with 
intramolecular vibrations is: 



(10) 
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4, respectively. Then, the value a — 0.28 is determined 
from the dipole-dipole coupling (-53 cm -1 ^) between 
sites 3 and 4. We constrain the reduced dimer model 
to match the derealization degree a — 0.28 and the 
experimentally determined exciton 1 and 2 energy dif- 
ference Sfafa — 158 cm -1 —, which results in renormal- 
ized site energies 8 = 135 cm -1 and coupling strength 
V = —41cm -1 . Given that our interest is the 1-2 cross 
peak of FMO whose beating signal is a manifestation 
of excitonic coherences, Figj5ja) shows the result of the 
numerical calculation of excitonic coherence. In the ab- 
sence of electron- vibration coupling, the amplitude of the 
coherence should be monotonically decreasing, since a 
dimer model will only have a single excitonic frequency in 
the one excitation manifold. The most prominent mani- 
festation of the presence of the vibronic structure is ob- 
served here, as a modulation of the excitonic coherence 
amplitude. In order to provide a figure of merit beyond a 
qualitative estimation for this modulation, an analytical 
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H = H bl + H b2 + J2 w<V**i I*) (*l + J2 u * l z > (*l + y d 3 > < 4 I + I 4 ) < 3 D 
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FIG. 5. Excitonic coherence P2,i(i) = Tr{|i/)i) (ip2\ p(t)}- In (a) full numerical calculation (continuous) and electronic coherence 
analytical result (dotted). In (b) are shown the Feynman diagrams that contribute to oscillating terms at frequencies concerning 
the excited electronic states on the response function. In (c) the result of the homogeneously broadened signal is presented 
for the the addition of both paths, at a waiting time t% = 20 fs. In (d)-(e), the waiting time resolved signal at points 
{uii,u>2} = {£ty,;L,<Aty, 2 } and {ui^ 1 + lu,lj^ 2 } respectively, is presented. Thick continuous line represent the real part of these 
oscillatory contributions. For comparison purposes, we show the absolute value of the rephasing signal (dotted) and the 
amplitude envelope from electronic coherence solution eq. (|ll[l (thin, continuous) scaled in amplitude and shifted in phase as 
described in the text. 



study is much preferable. 

In the excitonic basis, the projector |3) (3| = 
^(IViXV-il + a 2 1^2)^21 +a|^)(V 2 | +a|V> 2 > 
with a similar expression for |4) (4| by the exchanges 
ipi — > —%p2 and "02 - ► ipi j shows that in situations 
where 5 is much greater than the Coulomb interaction, 
the Hamiltonian eq. ([TUl) is diagonal in the excitonic ba- 
sis. With this approximation the excitonic coherence 
Tr{|-0i) (-02 1 p(t)} can be calculated using the total den- 



,— iHt 
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sity matrix evolution p(t) — e p(0)e 
our interest lies in the dynamics on the waiting time do- 
main, we suppose an initial state p(0) which is a product 
of the ith oscillator thermal state p b i and an electronic 
superposition proportional to the site transition dipole 
moments, p e ~ V, , /',//,, \i) (j\ = Y,iHiN>i IV'i} to 
best describe the situation after the arrival of two broad 
band low intensity pulses within the single excitation sub- 
space. After tracing out the exciton manifold and then 
noting that the trace over the oscillators can be arranged 
as a product of traces due to commutation among oper- 
ators describing vibrations in different chromophores, we 
make use of the displacement operator Ui — exp[— Si(bi — 
&■)] property: U l (H bi + v /s7u; l X l )[/,f 1 = H b i - s l u> l . 
This procedure leads to the interaction picture displace- 
ment operator f7j(i) = e lHbit Jj \e~ %Hhit correlation func- 
tion Trbi {U^ 1 (— i)UiPbi } , which for harmonic oscillators, 



is a standard result. Hence follows: 

Tr{\1>i) (Hp(t)} - jr^| e -^ 2 t+ 93 (t)+ 9 :(t) 

5^, 1 ^ 2 t — 2s coth(/3cj/2)(l— cos(uit)) (11) 



where gi(t) are the linewidth functions (remind eq.®) 
addressing the properties of vibrations in the ith site. 
As shown in the second line of the above equation, if 
oscillators are identical, the addition of the linewidth 
functions exp[g3(£) + gl(t)] results in a real argument 
exp[— 2scoth(/3Sw/2)(l — cos(wf))] presented (dashed) in 
FigJSJa), that quantitatively reflect the amplitude mod- 
ulation of the excitonic coherence. Also the complete 
solution, eq. pT|) in dotted lines shows excellent agree- 
ment with full numerical simulations. This underline the 
fact that disregarding excitonic population transfer terms 
mediated by vibrational modes cx jT^\ipi) {ipj\ seems 
correct in the B825, for the description of the vibronic 
wave packet motion. Moreover, this agreement in the 
current situation where V ~ S, highlights that in a more 
general context, the periodic modulation of the excitonic 
coherence, i.e., of the signal amplitude, is a signature of 
coherent dynamics occurring in the quantized vibrational 
manifold. The analytical result shows that the modula- 
tion of amplitude in the excitonic coherence will present 
a frequency equal to that of the prominent vibrational 
mode, namely uj, and will have a depth proportional to 
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the coupling of such vibration with the electronic transi- 
tions, via the Huang-Rhys factor s. 

Therefore, for practical purposes, the assumption that 
the Hamiltonian eq. (fTU|) is diagonal in the exciton ba- 
sis when accounting on vibronic effects on the 1 and 2 
excitonic coherence in the FMO, is well supported. We 
proceed with this assumption to investigate whether the 
modulation of excitonic coherence amplitude can be read- 
ily observed in the 1-2 FMO cross-peak from 2D spectra. 
The diagrams in FigJ5]Jb) are the contributions for the 
excitonic cross peak that present excited state evolution 
in the waiting time at a frequency equal to the energy 
difference Cj^, 2 — uj 1 p 1 = S^ 1% /, 2 , with corrections arising 
from the line-width functions. These oscillating terms 



carry dynamical information that provide access to the 
interference of the most salient dynamical contributions, 
and read: 

S 12 (ti,t 2 ,t 3 ) = 2Im[exp[-ioty 1 *i + iu)^ 2 t 3 + i5^ 2 t 2 ] 
x (M^iM^2)"V2,/M/,V2 exp[/i(ii,t 2 ,*3)] 

+ ^l 1 l4 2 exp[f 2 (t 1 ,t 2 ,t 3 )})} (12) 

where n^ it f is the dipole moment of the transition \tpi) — ¥ 
|/), and f\{ 2 ){ti,t 2l t 3 ) correspond to functions associ- 
ated with the cumulant expansion of the linewidth func- 
tiond of the paths shown in Figj5{b) under the assump- 
tion of a diagonal exciton-phonon interaction^ 1 -: 



fi{ti,t 2 ,t 3 ) = -g^ 2 (t 2 ) - g}f{h) - g/^ih + t 2 ) + gl 2f {t 2 ) + g^ 2 f(t 3 ) + g^^ih +t 2 ) - sy^^i) - g% 2 ^{h) 

+9fi>i(ti + h + t 3 ) - ffViV-i (*i + t 2 + t 3 )+g^ l (t 2 + t 3 ) - g^ 2 f{h + t 3 ) + g* flPl (t 3 ), 
f 2 (h,t 2 ,t 3 ) = —g^^iti + t 2 ) — g^foih) + g^ 2 (t 2 ) + g^ 2 [ti +t 2 + t 3 ) - g^^fa) - .9^2 ih +h) 



These expressions involve photon-echo rephasing con- 
tributions with linewidth functions of the singly ex- 
cited exciton states \ipi) , IV^) and doubly excited states 
l e 3) e 4) = I/)- The transformation of the coordinates 
Xi in the exciton basis yields to Xf = X 3 + X^X^ = 
j^(X 3 + a 2 X 4 ),A> 2 = j^(a 2 X 3 + Xi), which al- 
lows calculation of the linewidth function g^ (t) in terms 
of the site linewidth function g(t). For identical vi- 
brations in each site, it can be written gij = djg(t), 
with coefficients {C//,C^/ = Cf^ , C^ i ^ lj , } = 

{2,l, IT g F , I ^ T } for i + j. When a -> 0, the 

above expressions fulfill fi{t\,t 2 t 3 ) — f 2 {t\,t 2 t 3 ), while 

/ViJ = i-H'i = ^4 anc ^ VH>3,f ~ A*V>i = ^3 1,0 produce the 
expected vanishment of non-diagonal peaks in absence of 
electronic coupling among sites. 

Following the same procedure used for the diagonal peak 
in the homogeneously broadened FT signal, we per- 
form a first order Taylor expansion of w and present 
in Figsj5jc)-(e) the result of the 2D non-diagonal peak 
following eq. (fT2|) . The 2D spectra at early waiting 
times in Fig(5jc), shows appreciable contribution from 
the side-band structure. At the coordinate {^1,1^3} = 
{uj^, ig ,ui t (, 2 g} (FigJSJd)) the modulation of the signal en- 
velope is absent and supports that its principal contri- 
bution will mainly arise from ZPL transitions. On the 
other hand, the signal arising from the satellite structure 
{u>i,uj 3 } = + u;,u;,/, 2 } shows an appreciable modula- 
tion of its amplitude. This modulation has a frequency 
~ ui, that can be accounted for in virtue of ea. ([TT1) to the 
vibrational wave packet motion. 

The modulation in the signal of the cross peaks makes 
it a suitable witness for discrete vibronic structure in the 
environment. We follow with the detailed study of the 
progression from homogeneous to inhomogeneous broad- 



ening of the signal, which leads to interesting conclusions. 



C. Ensemble inhomogeneities: On the static nature of 
the difference among site and excitonic coherence decay 

Two different scenarios are shown in FigJSJa) for the 
introduction of static disorder. In the first case (dashed 
line) , the signal is averaged assuming uncorrelated static 
fluctuations, distributed according to the product of 
probability density functions from individual energies 
P{Cji)P{Cj 2 ), using comparable inhomogeneous widths on 
each. In such a case, the damping of the oscillations oc- 
curs in a time-scale « 100 — 200 fs, which docs not repre- 
sent the experimentally measured resul t 22 i 42 . By virtue 
of the microscopic model proposed here, the homoge- 
neous decay presented in FigJS](d)-(e) is similar, while the 
decay due to uncorrelated variation in the inhomogeneous 
average in FigJBKa) dashed line, is far too fast, compared 
to the appreciable visibility experimentally found for os- 
cillations well beyond t 2 = 1 p a 22 ' 42 . Hence it follows 
that there must exist a degree of correlation in the static 
energy variations in order to reproduce the experimental 
evidenc o 22 i 42 . Such correlation has been recently outlined 
as a process of paramount importance for understanding 
the signal arising from macroscopic ensembles^. 
Studies from first principles with molecular dynamics 44 , 
suggest that the energy shift in individual BChls in the 
FMO are due to two contributions: 1) local, arising from 
different geometries of each BChl deformed within their 
fitting into the protein scaffold, 2) environmental, in- 
duced by the charge density of the protein due to its crys- 
tal structure. The latter process is generally accounted 
for by the charge density of a-helices which extend across 
spatially extended areas involving in general, neighboring 
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BChls. For example, BChla 3, 4 and 5, are situated close 
to the a-helix 5^. Therefore, the electric fields generated 
by the protein are created in spatially extended regions 
including several chromophores within a single monomer, 
and is expected that a set of neighboring chromophores 
will experience Stark shifts from conformational changes 
of a single a helix, i.e., a correlation on energy shifts. 
If the shifts share sign and approximate magnitude, this 
situation will lead within the ensemble to small deviation 
from the average in the difference uj^ 2 — uj^ 1 . 

The fact that the main contributions in excitonic states 
IV'i) and 1 1P2) come from the neighboring sites 3, and 4 
and 5 respectively, will imply in this scenario a degree of 
correlation in the energy of these excitonic states. Hence, 
FMO monomers with exciton energies uj^ )1 +e and u!u 2 +e' 
fulfilling (e' 2 ) ~ (e 2 ) s» a 2 E , will result in e ~ e' and 
((e — e') 2 ) -c cte- With this assumption, is straightfor- 
ward to average the signal according to the product of 
gaussian distributions P(e), P(e — e') which describe the 
correspondent variations of exciton 1 energy and of the 
relative difference among excitons 1 and 2 energies, hav- 
ing standard deviation a e — ue and <r c _ C ' <C a t . The 
result of the correlated averaging procedure is presented 
in Figj6ja) (continuous line) which shows an appreciable 
amplitude of oscillations in the signal at ti = 2 ps, in 
good agreement with experimental evidence. It is impor- 
tant to highlight that this model reproduces the inhomo- 
geneous width recorded in absorption spectra or in the 
wi and W3 dimension in 2D spectra, while concomitantly, 
it is able to reproduce the enhanced lifetime of oscillation 
in the waiting time domain. 

This last conclusion requires a more thorough analysis. 
Remember the Feynman diagrams FigJSJb), where the 
signal in the coherence time t\ evolves at a frequency cy- 
cling among ground and excited electronic states, while 
that in the waiting time promotes oscillations among dif- 
ferent excitonic ket-bra combinations. The local field 
will therefore lead in t\ domain to a short dephasing de- 
cay time (as happens with the dashed line in FigH^a), 
~ 100—200 fs) proportional to the magnitude of the Stark 
shift due to the static fluctuations. The correlation in the 
contributing chromophores of the excitons allows (as pre- 
sented in FigJBJa)) the existence of one order of magni- 
tude difference in the decay of the signal when compared 
to non-correlated inhomogeneity. Hence, regardless of 
the magnitude of the Stark shift but sensitive to the gra- 
dient of the field producing it, the dephasing rate in the 
waiting time ti will decay slower than in t\ when cor- 
related fluctuations in the Stark shifts are present. The 
extended local field induced Stark shift picture is fur- 
ther supported by the recent experimental finding^ that 
cross-peaks associated with excitonic states whose main 
contributions involve neighboring chromophores have sig- 
nals (in the waiting time domain) with smaller dephasing 
rates. It cannot be explained based on heuristic argu- 
ments that shifts should be correlated by having equal 
signs at a given time. A more detailed analysis based on 
the actual charge distribution over the a helices may ad- 



dress microscopically this possibility, but will not to be 
treated here. 

The oscillation of the signal when correlated static in- 
homogeneities are present in FigjBJa) , resemble extremely 
well the behavior of the signal in the FMO 1-2 cross peak 
reported i n 22 ' 42 . This has been achieved by doubling the 
Huang-Rhys factor of this representative vibration well 
within the bounds 0.12 < s < 0.5, s = 0.24, and decreas- 
ing the homogeneous width of the mode 7—^7/8 which is 
acceptable from the fact that fluorescence line narrowing 
or spectral hole burning already give an inhomogeneously 
broadened ensemble average of the vibrational manifolds. 
The lower bound s = 0.12 has been used to infer the 
minimal effect of the mode of interest and keep consis- 
tent with experimental evidence. Moreover, the discus- 
sion concerning the diagonal peaks enhancement due to 
inhomogeneity and phase shift remain valid for s = 0.24. 
Quantitatively, the enhanced Huang-Rhys factor results 
in an increased contrast of the vibronic bands compared 
to the ZPL. Experimentally, the modulation of electronic 
coherences has been ascribed solely to overlap from peaks 
1-2 and 1-3, which would produce a beating according to 
the recently estimated Hamiltonian^l, at a frequency of 
38 cm' 1 . However, this result highlights that such beat- 
ing could arise solely from the discrete contribution of the 
most prominent vibronic mode. Hence, it must be under- 
lined that many processes are observed at once in the 2D 
spectra and that beating of the amplitude in signals can 
arise from excitonic interference and electron- vibrational 
coupling. Recently, resolution of the waiting time dy- 
namics at different cross peaks has been achieved^ 2 -, and 
therefore, the data is available to explore whether an am- 
plitude modulation having a frequency sa 36 cm' 1 is a 
general feature. 

In order to clarify whether vibrational coherence is 
responsible for the long lasting oscillations present in 
the non-linear response of the FMO complex, structural 
modifications were implemented^ 2 -, whose effect was in- 
tended to change the vibrational manifold frequencies. 
Interestingly, none of the modifications produced appre- 
ciable changes in the electronic coherence dephasing rate, 
nor the oscillation frequency of the signal. Note how- 
ever that our result implies a change in the period of 
the modulating envelope, and not that of the underlying 
oscillations when structural modifications are generated. 
Moreover, the usual mechanism by which ensemble os- 
cillations arc damped from averaging dephased contribu- 
tions, does not induce dephasing among signals regarding 
modes of different frequency. Fig^b) shows that increas- 
ing the vibrational mode frequency inhomogeneities, will 
only decrease the beating depth without enhancing the 
dephasing rate of the signal. We model the signal averag- 
ing according to a gaussian probability density function, 
with (j^j of 10 cm' 1 in agreement withal. This figure 
shows that using a physically feasible Huang-Rhys fac- 
tor is enough to present an appreciable modulation of 
the signal amplitude whose depth decreases as the vi- 
brational mode inhomogeneities are larger. Accordingly, 



14 




FIG. 6. In (a), the result of averaging homogeneous contributions with independent energy variations equal to those reported^ 
of {<7(,,ov} = {102,80} cm -1 (dashed) and with excitonic energy variation with {a t , a e - e >} = {102,5} cm -1 (real part, 
continuous; absolute value dotted). In (b) the additional inhomogeneity arising from the vibronic ensemble is included, and 
the result of the averaging using gaussian distributions with standard deviations a w = {1, 10}cm _1 are presented in continuous 
and dashed lines, respectively. Dotted line is the absolute value for <r„ = 10cm -1 . In all plots s = 0.24 (see text). 



static vibrational mode heterogeneities will lead to a de- 
creased modulation of the signal amplitude but will not 
influence the overall decay of the signal. In this way the 
robustness to structural modifications of the signal can 
be understood due to this novel non-dephasing effect over 
the ensemble inhomogenities. 



IV. CONCLUSIONS AND PERSPECTIVES 

We proposed a theoretical model for the low energy 
band of the FMO complex comprising a vibronic discrete 
structure and study the possibility to detect such vibronic 
dynamics with non-linear 2D optical spectroscopy. For 
experimentally determined parameters that characterize 
the electron-vibration coupling in FMO chromophores, 
the model is consistent with linear and non-linear opti- 
cal spectra for the low energy band of the FMO com- 
plex. The main features of the vibrational side bands 
are variations of the signal phase in the neighborhood 
of diagonal peaks and a modulation of the amplitude in 
the waiting time resolved signal at cross peaks, with a 
frequency corresponding to that of the energy of the vi- 
brational discrete mode. Complementarily, the robust- 
ness of excitonic coherence to dephasing can be further 
enhanced by correlations in the static energy disorder 
within a single monomer, presumably related with spa- 
tially extended electric field inhomogeneities that lead to 
correlated Stark shifts involving nearest neighbor chro- 
mophores. 

We show that correlated static disorder can lead to 
sample inhomogeneities that induce about an order of 
magnitude enhancement of excitonic coherence lifetime 
as compared to the t\ signal decay. When correlated 
static disorder is present, the discrete vibronic structure 
is not sensitive to ensemble inhomogeneities. Only under 
this conditions, the behavior of coherences at the single 
molecule level will represent the ensemble averaged sig- 
nal. Lastly, the resilience of the waiting time domain os- 



cillations to structural modifications is understood in the 
1-2 FMO cross peak as vibrations generate a modulation 
in the signal amplitude, and not a variation in the signal 
intrinsic frequency. The signal presents a dichotomic sen- 
sitivity to excitonic and vibrational features, where the 
former are extremely sensitive to inhomogeneities while 
the latter are almost unaltered. Noticeable effects in the 
waiting time 2D signal are enhanced by both discreteness 
and coupling strength of the prominent vibrations. 

The fact that the mode here studied might resonantly 
couple other electronic transitions in the FMO, and 
therefore, be useful to enhance derealization among ex- 
citons, is a possibility that remains to be explored. Our 
study suggests that a detailed understanding of the vi- 
brational environment, even at low frequencies, is crucial 
to achieve a full understanding of the electronic dynamics 
of pigment-protein complexes. 
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